Enhanced Sampling Simulations

with Ensembler

In this notebook, we give examples for enhanced sampling simulations with Ensembler. We show how to execute, visualize, and analyze those simulations for both 1D and 2D systems.

The enhanced sampling technologies are briefly explained in order to directly use this notebook for teaching purpose.

Maintainers: [@SchroederB](https://https://github.com/SchroederB), [@linkerst](https://https://github.com/linkerst), [@dfhahn](https://https://github.com/dfhahn)

Loading Ensembler and necessary external packages

import os, sys
my_path = os.getcwd()+"/.."
sys.path.append(my_path)

import matplotlib
matplotlib.use('TkAgg')

import numpy as np
%matplotlib inline
from matplotlib import pyplot as plt

from ensembler.potentials.OneD import fourWellPotential, harmonicOscillatorPotential, addedPotentials, metadynamicsPotential
from ensembler.potentials.TwoD import harmonicOscillatorPotential as harmonicOscillator2D,  wavePotential as wavePotential2D, gaussPotential
from ensembler.potentials.TwoD import addedPotentials as addedPotentials2D, metadynamicsPotential as metadynamicsPotential2D
from ensembler.samplers.stochastic import langevinIntegrator, langevinVelocityIntegrator
from ensembler.system import system

##Visualisation
from ensembler.visualisation.plotSimulations import simulation_analysis_plot

Enhanced Sampling Simulations in 1D

Unbiased System

We start our walkthrough with an unbiased reference simulation of a four-well-potential and give a small introduction in how to set up a simulation with Ensembler. The four-well-potential is defined by the x-positions (\(a\)-\(d\)) and the y-position \(ah\) - \(dh\) of the four wells. If wished, the potential can be scaled in the y-direction using \(V_{max}\). Note, that the energy is given in units of \(k_BT\).

# Simulation Setup
pot = fourWellPotential(Vmax=4, a=1.5, b=4.0, c=7.0, d=9.0,  ah=2., bh=0., ch=0.5, dh=1.)

Here, we choose the stochastic Langevin integrator for the numeric integration of our system. The user has to set the step size \(dt\) and the friction coefficient \(gamma\). Note, that this integrator already contains a thermostat. The temperature of the simulation will be set during the system setup (see below).

# Simple Langevin integration simulation
integrator = langevinIntegrator(dt=0.1, gamma=15)

The system class wraps the integrator and the potential. Additionally, the initial position of the particle \(position\) as well as the temperature parameter \(temperature\) are set.

# Put Potential and Integrator together to generate the simulation system
system1 = system(potential=pot, sampler=integrator,  start_position=4.2,  temperature=1)

To start the simulation we define the number of steps and run sys.simulate(). The progress of the simulation is displayed by a progress bar.

#simulate
sim_steps = 2000
cur_state = system1.simulate(sim_steps, withdraw_traj=True, init_system=False)

print("Trajectory length: ",len(system1.trajectory))
print()
print("last_state: ", cur_state)
print(len(system1.trajectory))

Simulation:  Simulation:   0%|          | 0/2000 [00:00<?, ?it/s]

initializing Langevin old Positions

Simulation:  Simulation: 100%|██████████| 2000/2000 [00:04<00:00, 423.56it/s]

Trajectory length:  2001

last_state:  State(position=4.926048134912287, temperature=1, total_system_energy=3.3534870056942463, total_potential_energy=3.3534870056942463, total_kinetic_energy=nan, dhdpos=-7.163241787364573, velocity=None)
2001

After running the simulation, the simulation data can be displayed as a table using sys.trajectory. Note, that we used a position Langevin integrator that did not calculate the velocities explicitly. Therefore, the kinetic energy and velocity are not defined. If you want to calculate the velocities during the simulation use the langevinVelocityIntegrator instead.

system1.trajectory.head()

	position	temperature	total_system_energy	total_potential_energy	total_kinetic_energy	dhdpos	velocity
0	4.200000	1.0	0.160027	0.158622	0.001405	1.595959	-0.053017
1	4.066507	1.0	0.016497	0.016497	NaN	-1.595959	NaN
2	3.949941	1.0	0.008460	0.008460	NaN	-0.533124	NaN
3	3.732899	1.0	0.281340	0.281340	NaN	0.395094	NaN
4	3.708773	1.0	0.334720	0.334720	NaN	2.117321	NaN

The simulation results can be visualized using the build-in visualizing tool. The left panel displays the energy surface as well as all visited states. The start and end-state are colored in blue and red, respectively. The middle panel shows the probability density distribution of the simulation as well as a boxplot of the distribution. This plot can be used to check if the system was simulated successfully. The rightmost panel shows the development of the force over time.

Note that without enhanced sampling only the minimum around x=4 is sampled. It would require a very long simulation time to overcome the energetic barrier and hence a lot of computing time. Enhanced sampling methods were developed to speed up those slow processes. We will explore some of the most commonly used enhanced sampling methods in the subsequent notebook.

#plot
simulation_analysis_plot(system1, title="Position Langevin", limits_coordinate_space=list(range(0,10)))

(<Figure size 1152x288 with 3 Axes>, None)

Enhanced sampling

Enhanced sampling methods can be divided into time-independent and time-dependent methods. Time-independent biases stay the same throughout the whole simulation whereas for time-dependent enhanced sampling the bias is updated during the simulation time.

Time-independent bias

Umbrella sampling

Umbrella sampling is a time-independent enhanced sampling method. In this method energy barriers are overcome by adding a bias potential. Note that umbrella sampling requires the choice of a reaction coordinate along which the bias is added. The choice of a suitable reaction coordinate is non-trivial for high dimensional systems. For our low dimensional 1D we can simply chose the x-axis.

One of the most frequent umbrella sampling method uses hormonic potentials to restrain the sampling to a certain region of the potential. This is especially useful for sampling transition regions.

We start as in the unbiased case by defining the potential. However, we have to define two potentials; the original potential and the bias potential. The original potential is the same potential as defined above, the bias potential is a hormonic oscillator centered at \(x_{shift}\) and force constant \(k\). In this case we want to sample the transition region around \(x\) = 5.5. Therefore, we set the \(x_{shift}\) parameter to 5.5. The force constant \(k\) defines how much we constrain the system. The higher the energy barrier, the more constrain is needed.

To sample the full potential energy landscape, we can set up multiple simulations with different \(x_{shift}\) parameters. For subsequent analysis of umbrella sampling simulations (e.g. using the weighted histogram analysis method WHAM) it is important that the sampling space of the different simulations overlap. The higher the force constant \(k\), the more simulations are needed achieve the overlap.

#Simulation Setup
origpot = fourWellPotential(Vmax=4, a=1.5, b=4.0, c=7.0, d=9.0,  ah=2., bh=0., ch=0.5, dh=1.)
biaspot = harmonicOscillatorPotential(k=10, x_shift=5.5)

The addedPotentials function of the biasOneD class wraps any two 1D potentials together. Therefore, it is straightforward to generate the enhanced sampling system from the original and biased potential.

#Add the bias and the original system
totpot = addedPotentials(origpot, biaspot)

All subsequent steps are identical to the unbiased system. Note, that the starting position of the simulation has to match to the \(x_{shift}\) parameter (i.e. should be reasonable close to \(x_{shift}\)) in order to avoid starting the simulation at very high energy states.

integrator = langevinIntegrator(dt=0.1, gamma=15)
system2 = system(potential=totpot, sampler=integrator,  start_position=4.2,  temperature=1)

#simulate
sim_steps = 2000
cur_state = system2.simulate(sim_steps, withdraw_traj=True, init_system=False)

print("Trajectory length: ",len(system2.trajectory))
print()
print("last_state: ", cur_state)
print(len(system2.trajectory))
system2.trajectory.head()

Simulation:  Simulation:   0%|          | 0/2000 [00:00<?, ?it/s]

initializing Langevin old Positions

Simulation:  Simulation: 100%|██████████| 2000/2000 [00:04<00:00, 406.70it/s]

Trajectory length:  2001

last_state:  State(position=4.6790392050205085, temperature=1, total_system_energy=5.196657035979513, total_potential_energy=5.196657035979513, total_kinetic_energy=nan, dhdpos=2.327564924852771, velocity=None)
2001

	position	temperature	total_system_energy	total_potential_energy	total_kinetic_energy	dhdpos	velocity
0	4.200000	1.0	8.616871	8.608622	0.008249	-11.404041	0.128445
1	4.402588	1.0	6.666382	6.666382	NaN	11.404041	NaN
2	4.314500	1.0	7.420504	7.420504	NaN	7.772810	NaN
3	4.297875	1.0	7.578418	7.578418	NaN	9.349773	NaN
4	4.177137	1.0	8.874049	8.874049	NaN	9.647798	NaN

In order to visualize the resulting potential we again use the plotting function simulation_analysis_plot. The visualization shows that umbrella sampling samples the high energy transition region around x=5.5 very well. In contrast, the unbiased simulation (see above) was stuck at the minimum around x=4.

#plot
simulation_analysis_plot(system2, title="Position Langevin", limits_coordinate_space=list(range(0,10)), oneD_limits_potential_system_energy=[0,30])

[0, 30]

(<Figure size 602.656x150.664 with 3 Axes>, None)

Scaled potential

addedPotentials can add any potential classes with a given dimensionality. In the special case that the added potential is identical to the original potential but a prefactor, one can scale the potential. In the example below we chose a four well potential that we will scale down in order to cross the energy barriers more easily. The procedure to define the original and bias potential are the same as described above.

sim_steps = 2000

#Simulation Setup
origpot=fourWellPotential(Vmax=4, a=1.5, b=4.0, c=7.0, d=9.0,  ah=2., bh=0., ch=0.5, dh=1.)
# same potential but Vmax is different
biaspot = fourWellPotential(Vmax=-3.5, a=1.5, b=4.0, c=7.0, d=9.0,  ah=2., bh=0., ch=0.5, dh=1.)
#Add the bias and the original system
totpot = addedPotentials(origpot, biaspot)

integrator = langevinIntegrator(dt=0.1, gamma=15)

system3=system(potential=totpot, sampler=integrator,  start_position=4.2,  temperature=1)

#simulate
cur_state = system3.simulate(sim_steps, withdraw_traj=True, init_system=False)

print("Trajectory length: ",len(system3.trajectory))
print()
print("last_state: ", cur_state)
print(len(system3.trajectory))
system3.trajectory.head()

Simulation:  Simulation:   0%|          | 0/2000 [00:00<?, ?it/s]

initializing Langevin old Positions

Simulation:  Simulation: 100%|██████████| 2000/2000 [00:05<00:00, 378.05it/s]

Trajectory length:  2001

last_state:  State(position=4.236433910957732, temperature=1, total_system_energy=0.02775586294076469, total_potential_energy=0.02775586294076469, total_kinetic_energy=nan, dhdpos=-0.1961553085554619, velocity=None)
2001

	position	temperature	total_system_energy	total_potential_energy	total_kinetic_energy	dhdpos	velocity
0	4.200000	1.0	0.034941	0.019828	0.015114	0.199495	0.17386
1	4.138912	1.0	0.009497	0.009497	NaN	-0.199495	NaN
2	4.119039	1.0	0.006937	0.006937	NaN	-0.138721	NaN
3	4.231237	1.0	0.026544	0.026544	NaN	-0.118940	NaN
4	4.197230	1.0	0.019279	0.019279	NaN	-0.230543	NaN

#plot
simulation_analysis_plot(system3, title="position Langevin", limits_coordinate_space=list(range(0,12)))

(<Figure size 602.656x150.664 with 3 Axes>, None)

Through the visualization we get the confirmation that we scaled the original potential down to 12.5% of its original height (first panel, compare the y-axis). Accordingly, the system has now enough thermal energy to cross all energy barriers and all four minima can be sampled.

Time-dependent bias

Time-dependent biasing methods update the bias during the simulation time. A frequently used time-dependent method is metadynamics/local elevation. There, a gaussian potential is added to the positions that were already visited during the simulation. Therefore, visiting these positions again, is energetically less favorable then in the previous visit (energetic penal).

Note, that in case of a Gaussian bias potential, the mean of the Gaussian is set to the current position of the particle and its width should be chosen small enough to avoid a big penal for neighboring states. Step by step the energetic minima are “flooded” and the particle can cross barriers more easily. In most applications the bias is only added every \(n^{th}\) step to iterate between free diffusion and biasing steps.

Metadynamics / Local elevation

We first define the original four-well potential. To apply metadynamics/local elevation we use the metadynamicsPotential function. In the initialization we have to define the height (\(amplitude\)) and width (\(sigma\)) of the gaussian bias function added. This bias potential is added every \(n\_trigger\) steps to the current position.

Adding more and more potentials every step leads to an energy function that demands more and more computation time every step. To avoid slowdown of the simulation the metadynamic bias is usually stored and calculated grid based. This allow a much faster simulation but comes at the cost of additional input parameter and small numerical errors in the bias force calculations. To initialize the grid, the user has to define the minimum x-position (\(bias\_grid\_min\)) and the maximum x-position (\(bias\_grid\_max\)) the grid should cover as well as the number of grid bins. Note, that no bias will be added to values below \(bias\_grid\_min\) or above \(bias\_grid\_max\).

In our example system we want to sample all four energy minima. Therefore, it is sufficient to set the grid between 0 and 10, which covers all four minima.

For the metadynamics/local elevation simulation we will reduce the simulation length by the factor of 10. This makes it easier to visually distinguish, and thus understand, the effect of metadynamics/local elevation.

sim_steps = 200  # reduced simulation length

#Simulation Setup
origpot = fourWellPotential(Vmax=4, a=1.5, b=4.0, c=7.0, d=9.0,  ah=2., bh=0., ch=0.5, dh=1.)

#Performe metadynamics
totpot = metadynamicsPotential(origpot, amplitude=.3, sigma=.21, n_trigger=10,
                               bias_grid_min=0, bias_grid_max=10, numbins=1000)

integrator = langevinIntegrator(dt=0.1, gamma=15)

system4=system(potential=totpot, sampler=integrator,  start_position=4,  temperature=1)

#simulate
cur_state = system4.simulate(sim_steps, withdraw_traj=True, init_system=False)

print("Trajectory length: ",len(system4.trajectory))
print()
print("last_state: ", cur_state)
print(len(system4.trajectory))
system4.trajectory.head()

Simulation:  Simulation:   0%|          | 0/200 [00:00<?, ?it/s]

initializing Langevin old Positions

Simulation:  Simulation: 100%|██████████| 200/200 [00:01<00:00, 112.56it/s]

Trajectory length:  201

last_state:  State(position=3.6584160417432634, temperature=1, total_system_energy=1.9065232003899293, total_potential_energy=1.9065232003899293, total_kinetic_energy=nan, dhdpos=1.8767775504426218, velocity=None)
201

	position	temperature	total_system_energy	total_potential_energy	total_kinetic_energy	dhdpos	velocity
0	4	1	-0.001201	-0.001344	0.000143	0.0034275749807509046	-0.016898
1	4.01858	1	0.299954	0.299954	NaN	-0.00342757	NaN
2	4.02243	1	0.300597	0.300597	NaN	-0.107778	NaN
3	3.82179	1	0.320542	0.320542	NaN	-0.13845	NaN
4	3.8931	1	0.296069	0.296069	NaN	0.552529	NaN

Metadynamic systems are visualized with the default simulation_analysis_plot function. simulation_analysis_plot will display the resulting potential after the last simulation step.

#plot
simulation_analysis_plot(system4, title="position Langevin", limits_coordinate_space=list(range(0,10)))

(<Figure size 602.656x150.664 with 3 Axes>, None)

In the first panel of the visualization we can see how the system’s energy minimum around x=4 was flooded and the particle can cross neighboring energy barriers more easily. The longer one simulates, the flatter the whole energy surface become. Note however, that artifacts can arise once the system leaves the grid defined by \(bias\_grid\_min\) and \(bias\_grid\_max\). Therefore, these parameters have to be selected very carefully.

The Ensembler package also contains a build-in functionality to make short movies of the simulation using animation_trajectory. An example code is shown below.

%matplotlib inline
import tempfile
from IPython.display import HTML
from ensembler.visualisation.animationSimulation import animation_trajectory, animation_EDS_trajectory
#plot simulation
ani, out_path = animation_trajectory(system4, limits_coordinate_space=list([0,8]), title="Awesome MetaDynamics",resolution_of_analytic_potential=1000)

##put it into jupyter:

os.chdir(tempfile.gettempdir())
x = ani.to_jshtml()
HTML(x)

Once Loop Reflect

2D systems

Enhanced sampling can also be performed in 2D. In the subsequent examples we show for the same enhanced sampling methods as above, how to set up, bias and visualize 2D systems.

To keep the usage of Ensembler as simple as possible, using 2D systems only requires minor changes with regard to 1D systems. Therefore, the methods description of the 2D systems are kept short. A more detailed description can be found at the corresponding 1D system.

Unbiased system

We first simulate the 2D system without any bias. In order to perform a 2D simulation, one simply has to choose a 2D potential. The Ensembler package then automatically adjusts all other parameters. If one wishes to set initial coordinates or velocities for the simulation, those should be 2D array.

sim_steps = 2000

#Simulation Setup
origpot = wavePotential2D(amplitude=(100,100), multiplicity=[1/5.,1/5.], radians=True)

integrator = langevinIntegrator(dt=0.1, gamma=15)

system5=system(potential=origpot, sampler=integrator,  start_position=np.array([50,50]),  temperature=3)
system.conditions = []
#simulate
cur_state = system5.simulate(sim_steps, withdraw_traj=True, init_system=False)

print("Trajectory length: ",len(system5.trajectory))
print()
print("last_state: ", cur_state)
print(len(system5.trajectory))
system5.trajectory.head()

Simulation:  Simulation:   0%|          | 0/2000 [00:00<?, ?it/s]

initializing Langevin old Positions

Simulation:  Simulation: 100%|██████████| 2000/2000 [00:13<00:00, 146.55it/s]

Trajectory length:  2001

last_state:  State(position=array([47.35998922, 48.74101918]), temperature=3, total_system_energy=array(-194.70375326), total_potential_energy=array(-194.70375326), total_kinetic_energy=nan, dhdpos=array([-0.84446165, -5.92401848]), velocity=None)
2001

	position	temperature	total_system_energy	total_potential_energy	total_kinetic_energy	dhdpos	velocity
0	[50, 50]	3	-167.779	-167.81430581529048	0.035009	[10.880422217787398, 10.880422217787398]	[0.20779926867318976, 0.16382141844501952]
1	[49.70279414612564, 49.87194018146838]	3	-172.26361178659465	-172.26361178659465	NaN	[-10.880422217787398, -10.880422217787398]	None
2	[49.92023196273755, 49.788797720020746]	3	-170.8939431107266	-170.8939431107266	NaN	[-9.864285658764715, -10.44709538053651]	None
3	[50.001878533812075, 49.52048660972045]	3	-172.61760995830818	-172.61760995830818	NaN	[-10.611324618356527, -10.162072428726702]	None
4	[49.837094272019854, 49.055893733438566]	3	-178.26193591217236	-178.26193591217236	NaN	[-10.886726346675461, -9.223506908835738]	None

After the simulation we can visualize the simulation on the 2D energy surface. The left panel shows the energy surface using a heatmap and the sampled regions. The middle panel shows the corresponding sampling histogram. The right panels show the forces in x and y direction.

In the visualization, we can see that the system contains multiple energy minima and maxima. The simulation was started in a minimum. During the simulation time the simulated particle is unable to leave this start minimum.

simulation_analysis_plot(system=system5, limits_coordinate_space=[[0,90],[0,90]])

(<Figure size 1152x288 with 5 Axes>, None)

Umbrella Sampling

Umbrella sampling is performed with the addedPotential2D function that works analogous to the addedPotential function for 1D. However, it expects two 2D potentials as input. Here we make umbrella sampling with a harmonic potential centered around point [10,5]. If needed, one can define different spring constants in x and y direction.

sim_steps = 2000

#Simulation Setup
origpot = wavePotential2D(amplitude=(100,100), multiplicity=[1/5.,1/5.], radians=True)
biaspot = harmonicOscillator2D(k=np.array([10,10]), r_shift=np.array([10,5]))

#Add the bias and the original system
totpot = addedPotentials2D(origpot, biaspot)

integrator = langevinIntegrator(dt=0.1, gamma=15)

system6=system(potential=totpot, sampler=integrator,  start_position=np.array([50,50]),  temperature=3)

#simulate
cur_state = system6.simulate(sim_steps, withdraw_traj=True, init_system=False)

print("Trajectory length: ",len(system6.trajectory))
print()
print("last_state: ", cur_state)
print(len(system6.trajectory))
system6.trajectory.head()

Simulation:  Simulation:   0%|          | 0/2000 [00:00<?, ?it/s]

nDimensions
initializing Langevin old Positions

Simulation:  Simulation: 100%|██████████| 2000/2000 [00:15<00:00, 126.72it/s]

Trajectory length:  2001

last_state:  State(position=array([11.11980939,  6.64837039]), temperature=3, total_system_energy=array(-17.03559536), total_potential_energy=array(-17.03559536), total_kinetic_energy=nan, dhdpos=array([3.59708069, 1.80317402]), velocity=None)
2001

	position	temperature	total_system_energy	total_potential_energy	total_kinetic_energy	dhdpos	velocity
0	[50, 50]	3	17957.3	17957.18569418471	0.106701	[410.8804222177874, 460.8804222177874]	[-0.43009617629068647, 0.1685773196386405]
1	[47.84258715024018, 47.27379380139479]	3	15896.7514362986	15896.7514362986	NaN	[-410.8804222177874, -460.8804222177874]	None
2	[45.47393574586834, 44.596997347830246]	3	13949.507919749325	13949.507919749325	NaN	[-381.29077175545814, -423.33746418087026]	None
3	[43.117072820646875, 42.35795506052497]	3	12334.305007795745	12334.305007795745	NaN	[-348.2586707434501, -386.28720044826366]	None
4	[40.92730524633992, 39.864563961185304]	3	10815.770825360642	10815.770825360642	NaN	[-316.8046219772066, -357.2742473049512]	None

In the visualization we see that the energy minimum of the resulting (original + biased) energy surface is the chosen point [10,10]. The system quickly leaves its start region, moves toward the center point, and samples this region intensively. Therefore, using umbrella sampling, the system was able to escape its starting minimum.

simulation_analysis_plot(system=system6, limits_coordinate_space=[[0,90],[0,90]])

(<Figure size 1152x288 with 5 Axes>, None)

Scaled potential

Just like in the 1D case, the addedPotential2D function can be used to scale potentials. Here we scale our system to lower the energy barriers between the minima.

sim_steps = 2000

#Simulation Setup
origpot = wavePotential2D(amplitude=(100,100), multiplicity=[1/5.,1/5.], radians=True)
biaspot = wavePotential2D(amplitude=(-99.5,-99.5), multiplicity=[1/5.,1/5.], radians=True)

#Add the bias and the original system
totpot = addedPotentials2D(origpot, biaspot)

integrator = langevinIntegrator(dt=0.1, gamma=15)

system7=system(potential=totpot, sampler=integrator,  start_position=np.array([50,50]),  temperature=3)

#simulate
cur_state = system7.simulate(sim_steps, withdraw_traj=True, init_system=False)

print("Trajectory length: ",len(system7.trajectory))
print()
print("last_state: ", cur_state)
print(len(system7.trajectory))
system7.trajectory.head()

Simulation:  Simulation:   0%|          | 0/2000 [00:00<?, ?it/s]

nDimensions
initializing Langevin old Positions

Simulation:  Simulation: 100%|██████████| 2000/2000 [00:13<00:00, 146.86it/s]

Trajectory length:  2001

last_state:  State(position=array([42.75113903, 47.01865092]), temperature=3, total_system_energy=array(-0.82056029), total_potential_energy=array(-0.82056029), total_kinetic_energy=nan, dhdpos=array([ 0.07310689, -0.00155753]), velocity=None)
2001

	position	temperature	total_system_energy	total_potential_energy	total_kinetic_energy	dhdpos	velocity
0	[50, 50]	3	-0.798273	-0.839071529076449	0.040799	[0.054402111088938554, 0.054402111088938554]	[0.16887982962108833, 0.23038403956280878]
1	[50.06819330012868, 50.475730212993724]	3	-0.8075835206686861	-0.8075835206686861	NaN	[-0.054402111088938554, -0.054402111088938554]	None
2	[50.00424684578778, 50.37561058471167]	3	-0.8172423122424846	-0.8172423122424846	NaN	[-0.055541397070397736, -0.06212744554758132]	None
3	[49.86712702854522, 50.34394534440884]	3	-0.8264623008381022	-0.8264623008381022	NaN	[-0.0544733596045468, -0.060546034704390905]	None
4	[49.447518272221835, 50.519819706148084]	3	-0.8360144240351701	-0.8360144240351701	NaN	[-0.05215336648225666, -0.0600407921407804]	None

simulation_analysis_plot(system=system7, limits_coordinate_space=[[0,90],[0,90]])

(<Figure size 1152x288 with 5 Axes>, None)

In the visualization we can confirm that the energy barriers are now low enough for the system to escape its initial minimum.

Metadynamics / Local elevation

Metadynamics/Local elevation changes the energy surface of the simulation during the simulation time (time dependent bias). Similar to its 1D analog, metadynamicsPotential2D uses a grid to store the bias added to the system. The minimal and maximal values of the grid in x and y direction have to be defined by the user. The gaussians added by the metadynamic bias can have elliptical shapes. The shape is defined by the \(sigma\) parameter. Again we use 10 times less steps then in the previous simulations.

sim_steps = 200

#Simulation Setup
origpot = wavePotential2D(amplitude=(100,100), multiplicity=[1/5.,1/5.], radians=True)

#Perform metadynamics
totpot = metadynamicsPotential2D(origpot, amplitude=10, sigma=(5,5), n_trigger=10, bias_grid_min=(0,0), bias_grid_max=(100,100), numbins=(1000,1000))

integrator = langevinIntegrator(dt=0.1, gamma=15)

system8=system(potential=totpot, sampler=integrator,  start_position=np.array([50,50]),  temperature=3)

#simulate
cur_state = system8.simulate(sim_steps, withdraw_traj=True, init_system=False)

print("Trajectory length: ",len(system8.trajectory))
print()
print("last_state: ", cur_state)
print(len(system8.trajectory))
system8.trajectory.head()

Simulation:  Simulation:   0%|          | 0/200 [00:00<?, ?it/s]

initializing Langevin old Positions

Simulation:  Simulation: 100%|██████████| 200/200 [00:05<00:00, 33.37it/s]

Trajectory length:  201

last_state:  State(position=array([46.07932953, 10.80526647]), temperature=3, total_system_energy=array(100.47096969), total_potential_energy=array(100.47096969), total_kinetic_energy=nan, dhdpos=array([-11.9848703, -17.3834667]), velocity=None)
201

	position	temperature	total_system_energy	total_potential_energy	total_kinetic_energy	dhdpos	velocity
0	[50, 50]	3	-167.748	-167.81430581529048	0.066539	[10.880422217787398, 10.880422217787398]	[-0.02531093277938583, -0.363919021799544]
1	[50.25308383817051, 49.48062679226069]	3	-70.1448261823999	-70.1448261823999	NaN	[-10.880422217787398, -10.880422217787398]	None
2	[50.51275537606412, 49.188256581477226]	3	-69.88844086011031	-69.88844086011031	NaN	[-11.727881288047744, -9.20424826668978]	None
3	[50.53329032897361, 49.129724189866835]	3	-70.0936788679769	-70.0936788679769	NaN	[-11.358214760746911, -9.342152134251553]	None
4	[50.493470641063915, 49.462597728123626]	3	-67.4766559535926	-67.4766559535926	NaN	[-11.422093346968953, -9.127151113762208]	None

From the visualization it is visible that regions which were already visited by the system show a higher energy. Especially the starting minimum is flooded with energy. Using metadynamics the system can rapidly sample the energy surface.

simulation_analysis_plot(system=system8, limits_coordinate_space=[[0,90],[0,90]])

(<Figure size 1152x288 with 5 Axes>, None)

This is the end of our Enhanced_Sampling Example notebook. We hope the reader has now a clearer understanding how Ensembler works and how different enhanced sampling methods can be used to increase the sampling.

Feedback is appreciated in our gitlab repository or directly to [@linkerst](https://https://github.com/linkerst).